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We investigate the dynamics of tlie Papapetrou equations in Kerr spacetime. These equations 
provide a model for the motion of a relativistic spinning test particle orbiting a rotating (Kerr) blaclc 
hole. We perform a thorough parameter space search for signs of chaotic dynamics by calculating the 
Lyapunov exponents for a large variety of initial conditions. We find that the Papapetrou equations 
admit many chaotic solutions, with the strongest chaos occurring in the case of eccentric orbits with 
pericenters close to the limit of stability against plunge into a maximally spinning Kerr black hole. 
Despite the presence of these chaotic solutions, we show that physically realistic solutions to the 
Papapetrou equations are not chaotic; in all cases, the chaotic solutions either do not correspond to 
realistic astrophysical systems, or involve a breakdown of the test-particle approximation leading to 
the Papapetrou equations (or both). As a result, the gravitational radiation from bodies spiraling 
into much more massive black holes (as detectable, for example, by LISA, the Laser Interferometer 
Space Antenna) should not exhibit any signs of chaos. 



PACS numbers: 04.70.Bw, 04.80.Nn, 95.10.Fh 



I. INTRODUCTION 

The present paper continues the investigation initi- 
ated in T|, which considered the dynamics of spinning 
test particles (i.e., compact astrophysical objects) orbit- 
ing a spinning black hole (Kerr spacetime). The primary 
objective is to determine whether the orbits of spinning 
compact objects spiraling into much more massive black 
holes are chaotic or not. Evidence for chaotic orbits in 
relativistic binaries has been presented in 0, 0, S IE IS 
(although only |3| was directly concerned with the ex- 
treme mass ratio systems we consider here). An astro- 
physical example of the systems we consider is a solar- 
mass black hole or neutron star orbiting a supermassive 
black hole in a galactic nucleus. Such systems are po- 
tentially important sources of gravitational radiation in 
frequency bands detectable by space-based gravitational 
wave detectors such as the proposed Laser Interferome- 
ter Space Antenna (LISA) mission. The methods of data 
analysis for signals from such detectors typically rely on 
matched filters, which use a discrete mesh in parame- 
ter space to achieve a given signal-to-noise ratio. The 
presence of chaos would cause the number of templates 
needed to fit a given signal to grow exponentially with 
decreasing mesh size, potentially overwhelming the com- 
puters performing the analysis. 

Many studies have used the Papapetrou equations to 
investigate the dynamics of spinning test particles in the 
background spacetime of a central black hole (including 
most recently [H S 0, H 113 ; see l9f and lldl for a more 
thorough literature review). We found in |2| that the 
Papapetrou equations formally admit chaotic solutions 
in Kerr spacetime, but the chaos seemed to disappear 
for physically realistic spins. This conclusion was ten- 



tative, since we investigated only a few values of the 
many parameters appearing in the equations. In the 
present study, we undertake a thorough search of param- 
eter space in order to determine the prevalence of chaotic 
orbits, both for dynamically interesting (but physically 
unrealistic) orbits with high spin parameter S and for 
astrophysically relevant systems with smaller spin. 

As in 0, we use Lyapunov exponents to measure the 
divergence of nearby trajectories and to provide an esti- 
mate for the timescale of the divergence. For chaotic sys- 
tems, initial conditions separated by a small distance eg 
diverge exponentially, with the separation growing as 
e(T) = £96^'': where A is the Lyapunov exponent. For 
nonchaotic orbits, A = 0, but when A is nonzero it pro- 
vides the e-folding timescale t\ = 1/ A for chaotic behav- 
ior to manifest itself. In this paper we use two rigorous 
methods for determining the maximum Lyapunov expo- 



nent Amax described in [J. Further details appear below 
in Sec. El and in Sec. Ill of Q. 

Our parameter space search makes use of a convenient 
orbital parameterization technique, which allows us to 
specify desired values of pericenter Vp and orbital inclina- 
tion i (Sec. IIII|I . We vary these orbital parameters along 
with spin magnitude, eccentricity, Kerr spin angular mo- 
mentum a, and spin inclination, and calculate the largest 
Lyapunov exponent for each set of initial conditions. A 
complete description of our methods for searching param- 
eter space appears in Sec.0 Though not exhaustive, the 
resulting survey of parameter space gives a thorough pic- 
ture of chaos in the Papapetrou equations with a Kerr 
background. 

We set G = c = 1, use sign conventions as in 
MTW ^-iid use Boyer-Lindquist coordinates {r,9,(f>) 
for Kerr spacetime. We use vector arrows for 4-vectors 
and boldface for Euclidean vectors. The symbol log de- 
notes the natural logarithm. 



'Electronic address: mhart l@tapir . caltech. edu] 
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II. SPINNING TEST PARTICLES 

We use the Papapetrou equations of motion 12] to 
model a spinning test particle. As reformulated by 
Dixon these equations describe the motion of a 

"pole-dipole particle," which neglects effects smaller than 
those due to the mass monopole and spin dipole (thus ne- 
glecting the tidal coupling, which is a mass quadrupole 
effect). The Papapetrou equations describe a particle 
with negligible mass compared to the masses responsible 
for generating the background spacetime. In our case, we 
write /i for the test particle's mass and M for the mass 
of the central Kerr black hole; the test particle limit then 
requires that /x <C M. 

A. Equations of motion 

Dixon's formulation for the equations of motion uses 
an antisymmetric spin tensor S^'^ to represent the spin 
angular momentum of the particle. The covariant deriva- 
tive of the 4-momentum deviates from the nonspin- 
ning (geodesic) case by a term representing a coupling of 
the spin to the background spacetime curvature. If we 
write for the 4-velocity, the full equations of motion 
appear as follows: 



^.p^^ = -iR^^^pv'^S'^P (2.1) 
VffS^" = 2p[^w'^]. 

Here R'^^^p is the Riemann curvature tensor of the space- 
time, which in our case corresponds to the Riemann ten- 
sor for the Kerr metric. Note that in the case of vanishing 
spin (all components of S'^'^ equal to zero) or flat space- 
time (all components of R^^j^^p equal to zero) , we recover 
the geodesic equation ^vP^ = 0. 

For numerical and conceptual purposes, we use a refor- 
mulation of the equations of motion in terms of the mo- 
mentum 1-form and spin 1-form S^. The spin 1-form 
can be derived from the 4-momentum and spin tensor 
using 

S^, = \e^,o.pu-S^P, (2.2) 

where u'^ — j [i. The spin 1-form is easier to visual- 
ize than the spin tensor, and the Papapetrou equations 
are better-behaved numerically in the 5' — > limit when 
expressed in terms of the spin 1-form. Details of this 
formulation appear in 1] and the appendix. 

B. General constraints 

The Papapetrou system of equations (|2.1|l is highly 
constrained, even in arbitrary spacetimes. The 4- 
momentum and spin satisfy normalization conditions 



that are conserved by the equations of motion: 

p% = -ii^ (2.3) 

and 

\Si^-S^,, = (2.4) 

A further condition is required to identify a unique world- 
line for the center of mass: 

V^.S^'' = 0. (2.5) 

A more detailed discussion of this "spin supplementary 
condition" appears in 0. 

C. Kerr constraints 

The Papapetrou equations share an important prop- 
erty with geodesic motion, namely, to each symmetry 
of the background spacetime (typically represented by a 
Killing vector) there corresponds a constant of the mo- 
tion. Kerr spacetime has two such symmetries, which 
provide two constraints in addition to those described in 
the previous section. (The Killing tensor that leads to 
the Carter constant in the geodesic does not correspond 
to a conserved quantity when the spin is nonzero; see 
Sec. IIIIB 31 below.) 

For each Killing vector ^ there corresponds a constant 
given by the standard expression for geodesies plus a con- 
tribution due to a coupling with the spin tensor: 

Q^C^m-Km;-^^"- (2-6) 

Kerr spacetime's two Killing vectors ^* = djdt and 
^■^ — d/d(f) then give energy and z angular momentum 
conservation: 

-E=pt~\9t^.,,S^''' (2.7) 

and 

Jz=P^-\g^^.,.S^'''■ (2.8) 

D. The spin parameter S 

The spin magnitude S that appears in Eq. H2.4|l quan- 
tifies the size of the spin and thus plays a crucial role 
in determining the behavior of spinning test particle sys- 
tems. As in we measure distances and times in terms 
of M and momenta in terms of /i. In these units, 5* is 
measured in units of /lAf . When measured in traditional 
geometrized units, the spin of a black hole can be as 
large as its mass squared, i.e., S'geom.max = M^- In this 
case the spin parameter goes like S ~ /x^/(/xM) = ^/M, 
which is necessarily small for a test particle system. This 
result applies to all astrophysically relevant systems, as 
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shown in 0]. Thus, we have the important physical con- 
straint that the spin parameter S must satisfy S <^ 1 for 
all physically realistic systems. A thorough discussion of 
various possibiHties (including neutron stars and white 
dwarfs) in showed that realistic values of S for LISA 
sources fall in the range 10~''-10~''. 

The smallness of the spin's effect is not dependent on 
our choice to measure S in terms of fiM. If instead we 
measure S in terms of /x^, then the equation of motion 
for 5''^'' becomes (rewriting V{t as D/dr) 



d{T/M) ~ / ' 



(2.9) 



where we have factored out from each variable its corre- 
sponding scale factor. This gives 



dr 



M 



(2.10) 



so that in these units the spin's effective magnitude is 
suppressed by a factor of fi/AI. Measuring S in terms of 
fiM effectively absorbs the unavoidable smallness of the 
spin effect into the spin parameter itself. 



III. PARAMETERIZING INITIAL CONDITIONS 

The many constraints on the equations of motion lead 
to considerable complexity in parameterizing the initial 
conditions. We summarize here the primary parameteri- 
zation method described in JJ , and then discuss in detail 
a method for parameterizing initial conditions using the 
geometrical properties of the orbit. 



S^; since the r-9 part of the Kerr metric is diagonal, we 
may then easily derive the 1-form components Sr and 5*61. 

Having specified values for seven of the variables, we 
must now use the five equations that relate them. Since 
we measure the particle's momentum in terms of its rest 
mass /i, the momentum normalization relation is 



(3.1) 



which we use to eliminate 
normalization condition 



We then solve the spin 



(3.2) 



to eliminate St/,. The spin-momentum orthogonality con- 
straint and the two integrals of the motion then give 
three equations in the three remaining unknowns pt, p^, 
and St'. 



= p^S^g'^'' 
E = -Pt + yt^^uS'^" 
Jz — P4> — ^g4>fi,iyS'^'^ 



(3.3) 
(3.4) 
(3.5) 



(The first of these equations is equivalent to the condi- 
tion p^S^" — 0, as discussed in the appendix.) 

We solve Eqs. (|3.3|) - (|3.5|l with a Newton-Raphson root 
finder, which works robustly given good initial guesses. 
The terms ^gt^^MS^" and ^g^^^^S'^" are typically small, 
even in the physically unrealistic case of S' '--^ 1, so 
and are good initial guesses for their corresponding 
momenta. We typically use Se as the initial guess for St. 
In all cases, we verify a posteriori that all constraints are 
satisfied to within machine precision. 



A. Energy and angular momentum 
parameterization 

It is easy to parameterize the Papapetrou equations di- 
rectly in terms of the momentum and spin components, 
but more physically relevant parameterizations are more 
difficult. The parameterization method discussed in 1] 
solves for the initial conditions using the integrals of the 
motion E and . This method also serves as the founda- 
tion for a more sophisticated parameterization in terms 
of orbital parameters fSec lIIIB)) . 

The energy and angular momentum parameterization 
method proceeds as follows. If we think of the system 
in terms of the spin 1-form S"^, we have twelve variables, 
four each for position, 4-momentum, and spin. Since 
Kerr spacetime is static and axially symmetric, without 
loss of generality we may set the initial time and axial 
angle to zero: t = (f> = 0. (We use the proper time r as 
our time parameter, as discussed in Sec. IIVI below.) We 
then provide the initial values for r, 9, and Pr, together 
with the constants E, J^, and S. Finally, we provide two 
components of the spin vector in an orthonormal basis. It 
is easiest to specify the radial and 9 components S^ and 



B. Orbital geometry parameterization 

While the method described above is sufficient for de- 
termining a set of valid initial conditions, parameterizing 
orbits by energy and angular momentum is not particu- 
larly intuitive. It is much more natural to think in terms 
of the orbital geometry, so we prefer to parameterize by 
the pericenter r^, the eccentricity e, and the orbital incli- 
nation angle l. Such parameters only have precise mean- 
ing for geodesic orbits, but are nevertheless still useful 
for the case of spinning test particles. 

1. Kerr geodesies 

The first step in parameterizing solutions of the 
Papapetrou-Dixon equations using orbital parameters is 
to solve the geodesic case. The traditional method for 
specifying a geodesic in terms of conserved quantities 
uses the energy E, the z angular momentum L^, and the 
Carter constant Q ^|. In order to use the orbital pa- 
rameters, we adopt a mapping from (rp, e, (.) to {E, L^, Q) 
based on unpublished notes supplied by Teviet Creighton 
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and Scott Hughes (and implemented in Hughes's Kerr 
geodesic integrator 14]). 

In order to use the more intuitive orbital parame- 
ters, we must determine the set {E,Lz,Q) given the set 
(rp,e,L). We obtain two of the necessary equations by 
noting that the radial velocity dr/dr vanishes at pericen- 
ter and apocenter, since these radii correspond to turning 
points in the radial motion. The equation for the time- 
evolution of the Boyer-Lindquist radius r is jlH : 



R(r) 



(3.6) 



where 

R{r) = [E{r^+a^)-aL,]^~A[r^ + {L,-aEf+Q], (3.7) 
and we use the standard auxiliary variables 



and 



cos^ I 



2Mr + a^ 



(3.8) 



(3.9) 



The quantity a is the Kerr spin parameter J/M, i.e., 
the central black hole's spin angular momentum per unit 
mass, which is dimensionless in our normalized units. 
From Eq. (|3.7|l we see that dr/dr = implies that 
R{r) — 0, so we obtain one equation at each turning 
point :^ 



R{rp) = 



and 



R{ra) - 0, 

where the apocenter is defined by 

1 + 



1 - e 



(3.10) 



(3.11) 



(3.12) 



The final equation required to complete the mapping 
is [3 



Q ^ l1 iav? L. 



(3.13) 



The value of l resulting from this definition agrees closely 
with the maximum value of \tt/2 — 6\ for a numerically 
determined solution to the equations of motion, i.e., it 
faithfully captures a geometric property of the orbit. 

Eqs. H3.10|I ~ H3.13() give three equations in three un- 
knowns, which are easy to solve using a nonlinear root 



^ In this paper we never consider exactly circular orbits, but we 
note that our prescription fails in this case: the conditions I3.1UI 



and 13. Hi are identical when e = 0, since Va 



For exactly 



circular orbits one must use the additional condition R'{rp) = 0, 
where R' = dR/dr. 



finder as long as good initial guesses for the energy, angu- 
lar momentum, and Carter constant can be found. The 
approach we adopt uses the degenerate cases of circular 
equatorial orbits to provide the raw material for analyt- 
ical guesses. The energies for prograde and retrograde 
circular orbits in the equatorial plane are 



i;P™(r) = 



Vl - 3w2 + 25i;3 



and 



1 - - dv^ 
Vl - 3t)2 - 2l^' 



(3.14) 



(3.15) 



where we write v = ^ Mjr and d — a/M for notational 
simplicity. The initial guess for the energy is then an 
average of these energies, weighted using the inclination 
angle, ^ with "radius" given by the semimajor axis of an 
ellipse with pericenter Tp and eccentricity e: 

= 1 [a+ £;P™(r,e„.i) + «- i?'-°'(^scmi)] , (3.16) 



where 



and 



a± = 1 ± cos t, 



1 - e 



(3.17) 



(3.18) 



It is possible (though rare) for Eq. H3.16|l to yield an 
energy guess greater than 1; in this case, we simply set 

Once we have a guess for the energy, we can find the 
corresponding guess for the angular momentum. Using 
the value from Eq. (|3.16l) and the expression for the an- 
gular momentum for a circular equatorial orbit gives 



L?"*'*''' = cos/ 



l-e2 



2(1 - £;gu°s>5-) 



as an initial guess for the angular momentum, 
the guess for the Carter constant is 



{LI 



f tan^t. 



(3.19) 



Finally, 



(3.20) 



Plugging Eqs. (|3.16() . (|3.19() . and H3.20|l into the nonlinear 
root finder yields the actual values of i?, L^, and Q to 
within machine precision in fewer than 10 iterations. 

One caveat about our parameterization method is 
worth mentioning: some values of (rp,e, i) correspond 
to unstable Kerr orbits, and in this case the method re- 
turns a set (E,Lz,Q) corresponding to an orbit with a 



^ We adopt the convention that a > 0, so that t indicates whether 
the orbit is prograde (0 < t < ^) or retrograde < t < tt). 
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different pericenter from the one requested. We can il- 
lustrate this behavior by factoring Eq. which is a 
quartic function in r: 

R{r) = (1 - E^){n - r){r - r2)(r - r3)(r - r^). (3.21) 

The roots are ordered so that ri > r2 > > r^. Bound 
motion occurs for ri > r > r2, which implies that ri = 
ra and r2 = Tp, but this only works for stable orbits. 
In the case that the orbit requested is unstable, the set 
(rp,e,t) returned by the algorithm instead corresponds 
to a nearby, stable orbit. In this case, the numerically 
calculated roots of R{r) satisfy rs = r^^ requested, but the 
parameterization method returns r2 as the pericenter, 
i.e., it returns a nearby stable orbit with pericenter = 
7*2 > ?'p, requested- As a lesult, the Qctual periccntcr is 
larger than the value requested. 

We use routines from to identify the boundary be- 
tween stable and unstable orbits (so that the latter may 
be excluded), but the code has a few minor bugs, and 
the identification procedure is not infallible. As a result, 
some such orbits appear in the results below (Sec. El 
and can be identified by having pericenters different from 
those requested. It is essential to understand, however, 
that the orbits returned by the parameterization algo- 
rithm are never unstable, and represent perfectly valid 
(stable) solutions to the equations of motion. 

2. From Kerr geodesies to Papapetrou initial conditions 

Once we have the set {E, L^, Q) for the geodesic case, 
we can use the well-known properties of the Kerr met- 
ric to solve for all the parameters necessary for the 
method described in Sec. IIII Al above. The first step is 
to force the constants of the motion to agree by setting 
-Bpapapetrou = -Ekoit and J^ = L^. Ncxt, wc set 6*0 = 7r/2 
(corresponding to the equatorial plane), since this is the 
only angle guaranteed to be shared by all Kerr geodesic 
orbits.'^ Finally, we must choose an initial value tq for 
the Boyer-Lindquist radius, which in the case of 9q = ti/2 
coincides with the crossing of the equatorial plane. One 
possibility is simply to use the average value of the peri- 
center and apocenter, 

ro^\{rp + ra). (3.22) 

This radius is guaranteed to lie between pericenter and 
apocenter, and the prescription in Eq. (|3.22l) works fine 
for the mildly eccentric orbits considered in this paper, 
but a more robust method must take into account that 



^ This restricts our sample to Papapetrou orbits that cross the 
equatorial plane. This almost certainly represents the vast ma- 
jority of valid Papapetrou solutions, but there remains the in- 
triguing possibility of spinning particle solutions that orbit per- 
manently above or below the equatorial plane. We leave an ex- 
amination of this possibility to future investigators. 



highly eccentric orbits should in general cross the equa- 
torial plane near pericenter. (Requiring a plane-crossing 
far from pericenter would force the inclination angle to 
be small, which is a constraint we do not wish to im- 
pose.) This suggests choosing an initial value of r close 
to pericenter. One fiexible method is to choose 

ro=rp(l + ae), (3.23) 

where a is a number of order unity. This reduces cor- 
rectly to Tp in the e = (circular) limit, and selects an 
equatorial plane crossing near pericenter in the e ~ 1 
limit. If we set a = 2 in Eq. |jn2Sll, then this method 
agrees exactly with Eq. H3.22|) in the cases e = and 
e = 0.5, and differs from Eq. by less than 15% 

when e = 0.6. Most of the results in this paper use 
Eq. H3.22|l . but Eq. H3.23|l is preferable in general. 

Once the initial 6 and r are known, we can determine 
all components of the Kerr 4-momentum p^, but three 
of the Papapetrou momenta are determined by the con- 
straints fSec. IIIlX)) . We are therefore free to force only 
one of the four components of the Papapetrou momen- 
tum pp to be the same as its Kerr counterpart. For the 
bulk of this paper, we set the radial components equal 
(p^ = Pp), since it is the radial momentum that is most 
closely tied to the stability and boundedness of the orbits. 
This choice results in Papapetrou orbits with pericenters 
and eccentricities fairly close to the Kerr geodesic values, 
but with much higher orbital inclination angles (Fig.EJ. 
The alternate choice of — pp results in Papapetrou or- 
bits with inclinations similar to their Kerr counterparts, 
but with very different pericenters (Fig. |SJ). See Sec. IVl 
for further discussion. 

We determine a value for the radial Papapetrou mo- 
mentum Pp using the Kerr geodesic parameters Ek , , 
and Q by applying Eq. (|3.7|l and the equation 

pV = (3.24) 

where p = r'^ + cos^ 6 • We then convert to Pr 
using Pr = p^grr- Proceeding exactly as in Sec. IIII Al 
we specify two of the spin components and eliminate two 
variables using 4-momentum and spin normalization, and 
then solve numerically for pt, p^, and St- The result is 
a set of initial conditions for the Papapetrou equations 
with the same energy and angular momentum as a Kerr 
geodesic with the desired values of Vp, e, and l. 

It sometimes happens that the Papapetrou initial con- 
ditions derived in this manner specify an orbit that is 
unstable against plunge into the black hole. Since there 
is no "effective potential" for a generic Papapetrou orbit 
as there is for Kerr geodesies, there is no way a priori 
to detect this instability. Plunge orbits are detected at 
runtime by testing for radial coordinates less than the 
horizon radius.* Orbits that plunge are removed from 
consideration since by definition they cannot be chaotic. 



* In practice, the most common runtime error is actually a nu- 
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3. Empirical orbital parameters 

In making the transition from geodesies to solutions of 
the Papapetrou equations, we are able to enforce the con- 
ditions i^papapctrou = £^Kcrr and = L^, but this is uo 
guarantee that the Papapetrou orbit has the correspond- 
ing orbital parameters rp and l: the spin-coupling term 
— ^^fiiuS'^" [cf. Eq. H2.6|l ] has a potentially large effect on 
the empirical values of the pericenter and orbital inclina- 
tion. This effect is most pronounced when we consider 
high spin parameter values, i.e., S ^ I. In these dy- 
namically interesting cases, the empirical pericenter and 
inclination will differ in general from the values requested 
in the parameterization. 

The empirical Papapetrou pericenter rp, p is easy to 
find: we simply integrate the orbit with a small step- 
size for a large number of periods, and then record the 
minimum radius achieved. In practice, this works ro- 
bustly, reproducing almost exactly the requested Kerr 
value of rp^ K in the limit <C 1. The only exception 
involves values of r^^ k that correspond to requested un- 
stable orbits, as discussed in Sec. IIII B II above. Each 
of these orbits has an empirical pericenter larger than 
the pericenter requested: the requested pericenter cor- 
responds to the root in Eq. H3.21|l . but the empirical 
pericenter returned by the algorithm corresponds to the 
larger root r2. 

Having found the empirical pericenter for an orbit, we 
must next find its empirical orbital inclination angle Lp. 
In order to reproduce the definition from Eq. I|3.13(l . we 
need to find an analogue of the Carter constant Q for 
spinning test particles. Kerr spacetime has a Killing ten- 
sor Kfj^i, 10 that satisfies 

= 0, (3.25) 

which gives rise to an extra conserved quantity in the 
case of geodesic motion: 

K = K^,p^'p\ (3.26) 

This quantity is related to the traditional Carter con- 
stant Q bv [nil 

K = Q+{L^-aEf. (3.27) 

When the test particle has nonzero spin, the quantity 
defined by Eq. H3.27|l is no longer constant, but there is 
an analogous expression that is conserved to linear order. 
Adapting a result from we can write this approxi- 
mately conserved quantity as 

C = K^.py' - 2p^S'>%n Up^ - fp fp..), (3.28) 
where 

U,^2acos9elel^+2relel^, (3.29) 



merical underflow in the integration stepsize as the particle ap- 
proaches the horizon. 



= 6 {^^elelel, + /fef.e^e^,^ , (3.30) 

and the {e°^} are the standard orthonormal tetrad for the 
Kerr metric. The effective Carter "constant" for spinning 
particles is then 

Q,s = C-{J,-aEf, (3.31) 

where we use the full angular momentum Jz (which in- 
cludes the contribution from the spin) in place of Lz ■ The 
quantity Qoff is nearly but not exactly not constant, so in 
order to define an empirical inclination angle we find the 
maximum effective Q over an orbit, and then define ip 

by 

Qeff.max = Jz t^U^ ^P- (3.32) 

As in the geodesic case, the value of t obtained from 
Eq. (|3.32|) agrees well with the maximum value of \Tr/2—0\ 
over an orbit. ^ When S = 0, Eq. reduces to 

the definition of the orbital inclination for geodesies, 

Eq. 12321). 



IV. LYAPUNOV EXPONENTS 
A. The principal exponent 

The Lyapunov exponents for a chaotic dynamical sys- 
tem quantify the chaos and give insight into its dynam- 
ics (revealing, for example, whether it is Hamiltonian or 
dissipative) . For a dynamical system of n degrees of free- 
dom, in general there are n Lyapunov exponents, which 
describe the time-evolution of an infinitesimal ball cen- 
tered on an initial condition. This initial ball evolves into 
an ellipsoid under the action of the Jacobian matrix of the 
system, and the Lyapunov exponents are related to the 
average stretching of the ellipsoid's principal axes. We 
described in jl] a general method for calculating all n of 
the system's exponents, but in the present study we are 
interested only in the presence or absence of chaos in the 
Papapetrou system, so we need only calculate the princi- 
pal Lyapunov exponent, i.e., the exponent corresponding 
to the direction of greatest stretching. A nonzero princi- 
pal exponent indicates the presence of chaos. 

When studying a differentiable dynamical system, we 
typically introduce a set of variables y = {yi} to rep- 



^ This is true only when we force the Kerr and Papapetrou values 
of p"" to agree ('Sec. llllB2l . which is the case for all orbits con- 
sidered except for the initial conditions shown in Fig. |S] In that 
case we revert to the simpler method of finding the maximum 
value of \it/2 — 9\ over several orbits. 
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resent the system's phase space, together with the au- 
tonomous set of differential equations 



(4.1) 



which determine the dynamics. Associated with each ini- 
tial condition is a solution (or flow). The principal Lya- 
punov exponent quantifies the local divergence of nearby 
initial conditions, so any fully rigorous method neces- 
sarily involves the local behavior of the system, i.e., its 
derivative. For a multidimensional system, this deriva- 
tive map is given by the Jacobian matrix. 



(Df) 



(4.2) 



It is the Jacobian map that determines the time- 
evolution of infinitesimally separated trajectories. If we 
consider a point y on the flow and a nearby point y -I- (5y, 
then we have 



f(y + Sy) = —(y + Sy) = — + — — , 

ar cIt cIt 

so that the separation Sy satisfies 



(4.3) 



^=f(y-(-<5y)-^=f(y + <Sy)-f(y). (4.4) 

CLT ClT 



Since 



f(y + <5y)-f(y) =Df -Jy + OdlcJylp), 



(4.5) 



we can write the time-evolution of the deviation vector 
as 



djSy) 
dr 



nf-Sy + Oi\\Sy\\ 



(4.6) 



If we identify the "infinitesimal" deviation Sy with an 
element ^ of the tangent space at y, we effectively take 
the limit as 6y —^ 0; the equation of motion for ^ is then 



dr 



Df 1. 



(4.7) 



This equation describes the time-evolution of the separa- 
tion between nearby initial conditions in a rigorous way. 
We will refer to ^ as a tangent vector, since formally is 
an element of the tangent space to the phase space. 

If the Jacobian Df in Eq. (|4.7I) were some constant 
matrix A, the solution would be the matrix exponential 



^(r) = exp(Ar) • |o 



(4.8) 



For long times, the solution would be dominated by the 
largest eigenvector of A, and would grow like 



IIWII 



(4.9) 



where Amax is the largest eigenvalue.^ (Here we have used 
ll^oll = !■) Turning things around, if we measured the 
time-evolution of ^, we could find an approximation for 
the exponent using 



Amax — 



log||l(r) 



(4.10) 



The general case follows by considering Jacobian ma- 
trices that are time-dependent. In this case, we are un- 
able to define a unique principal exponent valid for all 
times, but there still is a unique average exponent that 
describes the average stretching of the principal eigenvec- 
tor. Our method is to track the evolution of a tangent 
vector as it evolves into the principal axis of the ellipsoid. 
If we use Teir) = ||$(t)|| to denote the length of the 
longest principal ellipsoid axis, the principal Lyapunov 
exponent is then defined by 



log [reir)] 



(4.11) 



We use an infinite time limit in this formal definition, 
but of course in practice a numerical approach relies on 
a finite cutoff to obtain a numerical approximation. 

The Jacobian method for determining the largest Lya- 
punov exponent involves solving Eqs. H4.1|l and H4.7|) as a 
coupled system of differential equations in order to follow 
the time-evolution of a ball of initial conditions. One pos- 
sibility is then to use Eq. H4.11|l to estimate the system's 
largest exponent. A related technique, which provides 
more accurate exponents, is to sample log [''e(T)] at reg- 
ular time intervals, and then perform a least-squares fit 
on the simulation data. Since log [re (r)] = At, the slope 
of the resulting line then gives an estimate for the Lya- 
punov exponent. This is the method we implement in 
practice. 

A less rigorous but still useful technique, which we 
call the deviation vector approach, involves solving only 
Eq. H4.1|l . but for two initial conditions: yo and yo +(5yo- 
If the solutions to these initial conditions a time r later 
are y and y', respectively, then the approximate principal 
Lyapunov exponent is 



An 



log||'5y||/||'5yo|| 



(4.12) 



where Sy = y' — y. This approach has a serious draw- 
back: no matter how small the initial size of the devi- 
ation, eventually the method saturates as the difference 
between y and y' grows so large that it no longer sam- 
ples the local difference between two trajectories.^ On 



The most common choice for the norm |[ • || is the Euclidean norm, 
but we use a shghtly different norm in the case of the Papapetrou 
equations (Sec. llV^ belowV 

It is possible to re-scale the deviation once it reaches a certain 
size, but this method is error-prone since it can depend sensi- 
tively on the precise method of rescaling. The constrained na- 
ture of the Papapetrou equations also presents difficulties for 
rescaling, as discussed in Sec. II V Bl 
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the other hand, because we need only solve Eq. (|4.1|l and 
not Eq. H4.7|) , the deviation vector method is significantly 
faster than the Jacobian method (by a factor of approx- 
imately 5 for the system considered here). We therefore 
adopt the deviation vector method as our principal tool 
for broad surveys of parameter space. The method for 
handling the saturation problem is discussed in Sec. lIV (Jl 
A comparison of the Jacobian and deviation vector 
methods appears in Fig. ^ It is apparent that the two 
methods agree closely until the deviation vector approach 
reaches the saturation limit. 



B. The Papapetrou case 

The discussion in the previous section was of a general 
nature, applicable to virtually any dynamical system de- 
scribed by differential equations. Here we describe some 
of the details needed to apply the general methods to 
the Papapetrou equations. In particular, we must discuss 
further the key ideas of constrained deviation vectors and 
phase space metric. 

For the Papapetrou system, the phase space vector y 
has 12 components: 

y = (t,r,fJ.,(l),Pt,Pr,Pe,P,p,St,Sr,Sf,,S4,). (4.13) 

The tangent vector ^ has one component for each vari- 
able. The set of equations Eq. H4.1|l is simply the Papa- 
petrou equations written out in full: 



S^, = -p^(^R*y'S^v^p^Ss)+r"0^S^vP, 

where we have used the formulation in terms of the spin 
1-form described in the appendix. The second necessary 
equation, Eq. H4.7|l . requires the Jacobian matrix. 



dx^ dx^ dx^ 



dx" dp I, 
dp^ 



dx"" 



(4.14) 



V dx'' dp,, as"^ / 



whose explicit form appears in 0. 

It is important to mention that the tangent vector £ — 
or, equivalently, the deviation vector 5y — cannot have 
completely arbitrary components. On the contrary, the 
deviation must be chosen carefully, in order to ensure 
that, given a point y that satisfies the constraints from 
Sec. Ill Bl the deviated vector y + 5y also satisfies the 
constraints. Otherwise, the relation f (y -I- i5y) — f(y) ~ 
Df • 5y + 0{\\5y\\'^) is not satisfied. (This is the prin- 
cipal complication in implementing a rescaled version of 
the deviation vector method: the rescaled vector would 



violate the constraint.) In practice, we are able to find 
a valid deviation vector by applying the same techniques 
used to satisfy the constraints in the first place fSec. lIII|l : 
for details, see 

One final detail is the notion of metric: implicit in 
the definition of the Lyapunov exponent, Eq. H4.11I) . is a 
metric on phase space used to calculate the norm of the 
tangent vector |. We adopt a metric introduced in [l^ . 
which involves projecting the deviation vector onto the 
spacelike hypersurface defined by the zero angular mo- 
mentum observers (ZAMOs). The projection is effected 
using the projection tensor P'^ — S^,,+U^U,y, where is 
the ZAMO 4-velocity. The spatial and momentum vari- 
ables are then projected according to — > ; 



Pi = P^Pp-^ 'S'm ^ S'i. 



P^iSf,. After the 



projection, we calculate the Euclidean norm in the 3- 
dimensional hypersurface. We note that while this pre- 
scription is convenient, and reduces correctly in the non- 
relativistic limit, the magnitudes of the Lyapunov ex- 
ponents are similar for several other possible choices of 
metric 0. 



C. Chaos detector 

Since we are concerned with calculating Lyapunov ex- 
ponents for a large number of parameter values, we use 
the (unrescaled) deviation vector method because of its 



speed. As mentioned in Sec. IIV Al this method has the 
property of saturation (as illustrated in Fig. which is 
ordinarily a problem, but here we use it to our advantage 
as a sensitive detector of chaos. 

Our method for determining whether a particular ini- 
tial condition is chaotic is to consider a nearby initial 
condition separated by a small vector i5yo (with norm 
ll^yoll typically of order 10"'' or 10~^) and then integrate 
until the system reaches 90% of the saturation limit, de- 
fined as a separation Sy with unit norm. If we write 
Te = ll'^yll/ll'^yoll I the approximate Lyapunov exponent 
satisfies 



log [re (r)] = At, 



(4.15) 



so that A is the slope of the line log [?'e(r)] vs. r. Satura- 
tion occurs when re = l/||(5yo||, so that the integration 
ends when log [?'e(T)] = — log (0.9 ||5yo||). We record the 
value of log [7'e(T)] at regular time intervals (typically ev- 
ery time T — 100 M in our case), and upon reaching 90% 
saturation perform a least-squares fit on the simulation 
data to determine the exponent. 

We note that the cutoff value of 0.9 is somewhat ar- 
bitrary and is the result of numerical experimentation. 
When using the (unrescaled) deviation vector approach, 
most of the chaotic systems saturate — that is, plots of 
log [re(T)] vs. T flatten out (Fig. 0) — when the separa- 
tion is of order unity, corresponding to radial separations 
of order M, velocity separations of order 1, and angu- 
lar separations of order 1 radian. The 90% prescription 
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FIG. 1: Rigorous Jacobian method compared to unrescaled deviation vector method for an S' = 1 particle in maximal (a = 1) 
Kerr spacetime. The vertical axis is the natural logarithm of the largest principal axis ri of the phase space ellipsoid; the slope 
is the principal Lyapunov exponent, Amax « 5 x 10""^ A4~^ . The unrescaled deviation vector method started with a deviation 
of size 10~^, and it saturates at ~16. This corresponds to a growth of e^® « 9 x 10®, which means that the separation has 
grown to a size of order unity. In conventional units, this indicates radial separations of order GM/c? and velocity separations 
of order c. The norm is calculated using the projected norm described in Sec. If V 51 




20 40 60 80 100 
7- (10^ M) 



FIG. 2: A comparison of chaotic and nonchaotic initial con- 
ditions. The slopes of the lines are the approximate Lyapunov 
exponents. Each initial condition shares the same values of 
S = 0.1, e = 0.6, and rp — 1.21. They differ only in incli- 
nation angle: t = 31° (chaotic) and t = 28.5° (not chaotic). 
Their respective Lyapunov exponents are A = l.Ox 10~^ Af-^ 
and A = 3.0 X 10"^ M'K See SecEHfor details. 

ends the integration before the growth flattens out, so 
that the numerical estimate for the exponent is still rea- 
sonably accurate. 

For the large maps of parameter space, we typically 
integrate up to Tfinai = 10^ M or saturation, whichever 
comes first. We choose this maximum time mainly for 
practical reasons: it is the longest integration possible 
in a reasonable amount of time. [We integrate as deeply 
as 10^ M (Sec. IV Gp for individual orbits, but such long 
integrations are impractical for more than a handful of 
parameter combinations.] Dramatically longer integra- 
tions are also not particularly useful, since the timescale 
for gravitational radiation reaction is on the order of 
tgw = M^/fi :18|, which for the most relevant LISA 
sources is 10^-10^ M (i.e., /i - lO^^-lQ-^M). Search- 
ing for chaos in such systems on a timescale longer than 
~10^Af is probably pointless, since the radiation reac- 




T(10^ M) 

FIG. 3: A chaos mimic: log [^^(r)] vs. r for an S = orbit. 
The size of the initial deviation vector is eo = 3.3 x 10"^. The 
value of log [)"e(T)] periodically rises up to the saturation level 
(shown as a horizontal line at 14.82, since e^*"^"^ eo = 0.9). The 
system's spin satisfies 5 = 0, and is hence fully integrable, 
which implies no chaos. We detect such spurious chaos by 
demanding several saturation points in a row for a positive 
detection. The corresponding orbit appears in Fig. 1221 



tion would dominate the dynamics in this case. Finally, it 
appears that chaos, when present in the Papapetrou sys- 
tem, manifests itself on relatively short timescales (10^- 
10^ M), or else not at all. The onset of chaos is marked 
by a qualitative change from power law growth (which 
appears as logarithmic growth on our plots of log [re(T)] 
vs. r) to exponential growth (which is linear on the same 
plots) . An example of two similar initial conditions giving 
rise to qualitatively different dynamical behavior appears 
in Fig. 12 

One important refinement to the technique described 
above is to require several 90% saturation points in a 
row before declaring the orbit to be chaotic. This is 
necessary because some nonchaotic orbits have very high 
amplitudes on plots of log [re (t)] vs. t (Fig. OJ. This 
phenomenon occurs mainly for orbits with many periods 
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FIG. 4: The difference Sr between the Boyer-Lindquist radii 
of two nearby trajectories for a chaos mimic. The growth in 
the separation is substantial, but not exponential. The initial 
conditions are the same as in Fig.|3 
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T(IO^M) 
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FIG. 5: The difference 5r between the Boyer-Lindquist radii 
of two nearby trajectories for a chaotic orbit. The initial 
condition is taken from the inner region of Fig. |^ (rp, k ~ 
2.0, ck = 0.5, LK ~ 10°, S = 1). On this linear scale, the 
separation seems to burst unexpectedly, in contrast to the 
relatively smooth linear growth for the nonchaotic orbit shown 
in FigH 



in the deeply relativistic zone near the horizon. Such 
orbits may reach "90% saturation" briefly as part of their 
oscillation, but quickly return to separations far below 
our chaotic cutoff. We therefore adopt the criterion of 
three 90% saturation points in a row (with a time T — 
100 M sampling interval) as a robust practical test for 
chaos. See Sec. IV Hi for further discussion. 

Our confidence of this method's robustness derives 
from comparing the method above to the Jacobian 
method for the same initial conditions. Since the Ja- 
cobian method does not saturate, the agreement of the 
two methods indicates that our procedure provides an 
accurate detector for chaos (as in Fig. IT^ below'). 



D. Implementation notes 

We integrate the Papapetrou equations on a computer 
using Bulirsch-Stoer and Runge-Kutta integrators imple- 



mented in the C programming language, as described 
in 0. The derivatives and Jacobian matrix are exten- 
sively hand-optimized for speed. We monitor errors us- 
ing constraints and conserved quantities, with a global 
error goal of 10"^'^. The errors are at the 10~^^ level or 
better for highly chaotic orbits after 10^ M . Orbits with 
low spin or high pericenter are even more accurate, often 
achieving the error goal of 10" 

The many plots in Sec. El are typically generated using 
driver programs written in the Perl programming lan- 
guage, which in turns calls the C integrator repeatedly. 
This general paradigm — using an interpreted language 
such as Perl to call optimized routines in a compiled lan- 
guage such as C — is one we warmly recommend. 



V. RESULTS 

We present here the results of parameter variation in 
the Papapetrou system of equations. We represent the 
effects compactly using several different kinds of plots, 
most of which involve plotting inclination vs. pericenter, 
with other parameters held fixed. We refer to these as 
rp-L maps. As discussed in Sec. IIIII starting with the 
Kerr values rp ^i t/f, and ex, we find the correspond- 
ing energy Ek, angular momentum Lz^k, and Carter 
constant Qk, and then force the Papapetrou values to 
satisfy Ep — Ek, Jz,p = Lz^k and pp = p^. Each rp-L 
map has two components: part (a), shown on the left, 



uses the Kerr parameters lk and rp, k requested by the 
parameterization fSec. 1111 ti 2fl . while part (b), shown on 
the right, always shows the empirical Papapetrou values 
Lp and rp, p in the sense of Sec. 1111 B"3l 

One important feature of rp-i parameter space appar- 
ent in the empirical plots is the prevalence of large em- 
pirical inclination angles for all values of the Kerr incli- 
nation LK- Fig. 13 shows the mapping for = 1 between 
the requested Kerr parameters [Fig.[7|[a)] and the empir- 
ical Papapetrou parameters [Fig. jT^b)] using a shading 
scheme (which appears as a more informative color scale 
in electronic versions of this paper). We see that even 
the orbits at the bottom of Fig. |7{a) get mapped to high 
empirical inclination angles; lk — 1° orbits are mapped 
to inclinations of order ip = 40°. 

This compression of parameter space is the result of 
our choice to force the Kerr and Papapetrou value of the 
radial momentum to agree (Sec. 1111 B 2|) . The price we 
pay for this choice is that the 6 component of the Papa- 
petrou momentum — which is constrained to satisfy the 
equations in Sec. Ill Bl — is relatively large. This high pp 
flings even orbits with low requested values of lk to high 
inclinations (Fig. [7^. It is possible to find low-inclination 
Papapetrou orbits by requiring that the Kerr and Papa- 
petrou values of agree, but at the cost of forcing p"^ 
to be very different — again a result of the constraints. 
The resulting parameter space (Fig. ISJ is not nearly as 
compressed in inclination angle, but the Papapetrou peri- 
centers are compressed and shifted down, and many re- 
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FIG. 6: The orbit of a maximally spinning (S = 1) test particle in maximal (a = 1) Kerr spacetime, plotted in Boyer-Lindquist 
coordinates, (a) The orbit embedded in three-dimensional space, treating the Boyer-Lindquist coordinates as ordinary spherical 
polar coordinates; {h) y — r sin sin (f) ys. x — r sin 6 cos (j>\ (c) z vs. p = ya?^+ly2 . The requested orbital inclination angle is 
i = 6°, but the strong spin coupling gives rise to an empirical value closer to ip = 46°. The empirical pericenter is r-p = 2.219 M, 
which is fairly close to the requested value of 2.367 M. 



Geodesic (K-cn.") parameters 



Empirical Papapctioii paLametei-s 





FIG. 7: Shading/coloring of parameter space for an 5 = 1 orbit. The shading/color on the left is repeated on the right, so that 
we can visually determine the mapping of {rp,K,LK) to (rp^p,Lp). It is evident that orbits with low requested pericenters and 
orbital inclination angles are mapped to low-pericenter orbits at high inclinations, and the entire parameter space is compressed. 
Not that the gap in stable initial conditions visible in (b) is a true gap, not a fold. 




FIG. 8: Shading/coloring of parameter space for an S = 1 orbit illustrating the alternate parameterization method from 
Sec. nil B 21 As in Fig.|7| the shading/color on the left is repeated on the right. The spatial part of the initial spin is purely in 
the z direction. The Kerr and Papapetrou p% values are forced to agree, which leads to similar inclination angles in parts 
(a) and (b), at the cost of dramatically different pericenters. 
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quested values of r^, K axe lost as plunge orbits. In addi- 
tion, the empirical values of the eccentricity are typically 
not close to the value requested (reaching, e.g., ep — 0.75 
for ck = 0.5). Because of the deficiencies of this alter- 
nate parameterization method, we choose the fixed 
plots are our primary investigative tool in this paper. 



A. Varying pericenter and orbital inclination 

In this section we show r^-i maps for orbits with fixed 
eccentricity e = 0.5 and spin parameters = = 
0.2 5, where S is the total spin.^ We indicate the strength 
of the chaos at each point with a grayscale rectangle, 
with the darkest colors representing the strongest chaos 
(and with white indicating no chaos). An example of 
such a plot appears in Fig. O (Fig- El shows a simi- 
lar plot for the alternate parameterization method that 
forces = Pp. The maximum exponents in the two 
cases are virtually identical.) The most important gen- 
eral results evident from the plots are two-fold. First, 
the largest exponents occur for orbits with pericenters 
deep in the relativistic zone near the horizon. Second, 
the prevalence of chaotic orbits is a decreasing function 
of spin parameter S*, with virtually all chaos gone by the 
time 5" = 0.1. 

An example of the value of the empirical r^, p-Lp plots 
appears in Fig. 1261 which shows orbits of particles with 
spin S — 0.5. The appearance of a strongly chaotic point 
in Fig. 1261 seems perplexing, given that it is surrounded 
but many points with much smaller exponents. As is ev- 
ident from the empirical plot, this point of strong chaos 
(which is, in fact, the largest Lyapunov exponent for any 
of the plots) maps to a compressed area of initial condi- 
tions with small empirical pericenters [Fig. I26f b^]. 

From an astrophysical standpoint, the most interesting 
parameter to vary is the spin S, since only <C 1 orbits 
are physically realistic (Sec. Ill D|l . From Figs. I25H29I 
which involve varying S from 0.9 down to 10~*, we see 
that both the prevalence and strength of chaos decrease 
significantly as S is decreased. The Lyapunov exponent 
ranges as high as 10"^ M'^ when S = 0.5 (Fig.^B)), but 
the chaos is weak when S = 0.2 (Fig. I27|) and is gone 
below S = 0.1 (Figs. I2H1 and Onil. 



B. Varying eccentricity 

The choice of e = 0.5 in the previous section is par- 
tially motivated by likely gravitational wave sources for 
LISA 19], e.g., a neutron star or small black hole in an 
eccentric orbit around a supermassive black hole. In this 



These "fixed" spin parameters give rise to a variety of spin in- 
clination angles in the fiducial (ZAMO) rest frame; see Sec. EH 
below. 



section, we consider a second series of eccentric orbits at 
fixed e — 0.6 and varying spin parameter. We also inves- 
tigate the case of a near-circular orbit (e = 0.01) more ap- 
propriate for the circularized gravitational wave sources 
important for ground-based detectors such as LIGO. 

The rp-L plots for e — 0.6 follow the same pattern as 
those with e = 0.5. Chaos is strongest for orbits with 
small pericenters and values of S of order unity (Figs. 1201" 
I33|l . There is a single orbit at S* = 0.1 that appears to 
be chaotic (Fig. I32|l . but other than this one exception 
there is apparently no chaos below 5* = 0.1. A close 
examination of the single 5 = 0.1 chaotic orbit appears 
in Figs. I11H12I which shows that the chaos is real. 

The effect of chaos is smaller for the near-circular 
(e = 0.01) orbits considered fFigs. 01h37|l . with typical 
Lyapunov exponents of order 2 x lO^'^ M^^ when 5 = 1 
fFig l34|) . Moreover, we find only four points with nonzero 
exponents for S — 0.5 fFig. I35|) . in contrast to the more 
eccentric orbits, which have strong chaos for S — 0.5. 
By the time S — 0.1 (Fig. I36() . all chaos has completely 
disappeared for the near-circular orbits. 



C. Varying Kerr parameter a 

Here we investigate the effect of the Kerr parame- 
ter a on the strength and prevalence of chaos. Although 
Suzuki & Maeda found in that there is chaos even 
in Schwarzschild spacetime, to which Kerr spacetime re- 
duces when a = 0, the chaotic orbits found in 0] are 
exceptional orbits lying on the edge of a generalized ef- 
fective potential. We found evidence in fl| that such 
chaotic orbits are rare. 

The conclusion that chaotic orbits become less preva- 
lent as a —> is supported by a more thorough exami- 
nation, as illustrated in Figs. I38II43I All of these orbits 
have eccentricity e = 0.5 and spin S = 1, and a varies 
from 0.9 down to 0. Even for low values of a, the param- 
eter space is strongly affected by the spin, with plots of 
the empirical pericenter and orbital inclination showing 
significant distortion. Nevertheless, we see unambigu- 
ously that the chaotic orbits prevalent when a = 1 are 
greatly suppressed as a decreases, with no chaotic orbits 
at all below a — 0.2. This appears to be a result of the 
increase of the marginally stable radius r^s as a ^ 0. 
When = 0, the minimum stable radius is at r^g — 4 
in units of the central black hole's mass, fully 2M away 
from the horizon at rn = 2M. As discussed in the 
extra (spherical) symmetry of the Schwarzschild metric 
leads to an additional integral of the motion, which also 
has a suppressive effect on chaos. 



D. Spin cutoffs for chaos 

In this section we provide spin cutoff values for chaos, 
i.e., the minimum spin values for which chaos exists. For 
a given initial condition defined in terms of fixed orbital 
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FIG. 9: rp-L map for a = 1, e — 0.5, 5* = 1, 5"^ = 5*^ = 0.2 5": chaos strength as a function of pericenter and orbital inclination 
angle, (a) Requested parameters; (b) empirical parameters. The initial conditions are the same as in Fig. Q The rectangles 
are shaded according to the Lyapunov exponent for the initial condition represented by each point, with darker shades of gray 
representing larger exponents and hence stronger chaos. The point at r^, k ~ 2.3 and lk = 20° is one of the points in Fig. 1171 
which shows the effects of varying 5"^ and 5*^. The largest exponent in this plot is A = 4.1 x 10""^ , corresponding to a 
timescale of 1/A = 2.4 x 10^ M for a factor of e divergence in nearby initial conditions. 
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FIG. 10: Tp-i map for a — 1, e = 0.5, S — 1, using the alternate parameterization method that forces p%; — pp. The initial 
spatial component of the spin is purely in the z direction, (a) Requested parameters; (b) empirical parameters. The initial 
conditions are the same as in Fig. |H| The rectangles are shaded according to the Lyapunov exponent for the initial condition 
represented by each point, with darker shades of gray representing larger exponents and hence stronger chaos. The largest 
exponent in this plot is A = 4.0 x 10~^ , corresponding to a timescale of 1/A = 2.5 x 10^ M for a factor of e divergence in 
nearby initial conditions. These values are virtually identical to the values in Fig. |S] 



parameters (as described in Sec. 1111 Bll , we vary the spin 
parameter S and find the maximum value for which chaos 
occurs. The smaller this cutoff value, the stronger the 
chaos: nonchaotic orbits have a cutoff value of = 1, 
i.e., they are not chaotic even in the extreme S — 1 limit; 
conversely, a cutoff value of 5 = 10~^ would indicate 
chaos for the (physically realistic) value S = 10^^, but 
not for any smaller values.^ 



Fig. 1131 is an example of a spin cutoff map. The proce- 
dure for producing such a map is similar to the method 
used to make the Lyapunov Vp-L maps: we consider a grid 
of points in r^-t space, and for each point we find an ap- 
proximate value for the spin cutoff. We begin by finding 
out if the system is chaotic for S = 1, using the Lya- 
punov map as a start. If the orbit for S* = 1 is chaotic, 
we halve the spin value and calculate the Lyapunov ex- 
ponent for S = 0.5. If the system is still chaotic, we 
consider S — 0.25; otherwise, we consider S = 0.75. The 



Implicit in this scheme is an assumption of monotonicity, i.e., 
monotonically decreasing chaos as S decreases. While not strictly 
true (as discussed in Q), this assumption is still broadly appli- 



cable, and exceptions are rare. 
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FIG. 11: Natural logarithm of the largest ellipsoid axis vs. 
time. We use the unrescaled deviation vector method to inves- 
tigate the single chaotic initial condition from Fig. 1321 which 
has Tp, K — 1.21, ek = 0.6, lk = 31°, and S = 0.1. The prin- 
cipal ellipsoid axis grows until it saturates at t = 17600 M. 
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FIG. 12: Natural logarithm of the largest ellipsoid axis vs. 
time using the unrescaled deviation vector method and the 
rigorous Jacobian method. The unrescaled integration is iden- 
tical to that shown in Fig. 1111 The two methods agree un- 
til the deviation method saturates, at which point we stop 
the deviated vector integration. The continued growth of the 
Jacobian method confirms that the orbit determined by the 
initial condition is indeed chaotic. 



procedure continues until the difference between chaotic 
and nonchaotic spin values is smaller than some thresh- 
old, which we choose to be 0.05. (This value is chosen 
to achieve reasonable accuracy while still completing the 
calculations in a tolerable amount of time.) 

Figs. 1 131 1 151 show the spin cutoff values for several pa- 
rameter combinations corresponding to Lyapunov maps 
from Sec. IV Al above. The plots are color-coded with 
grayscale so that the most chaotic points — those with the 
smallest spin cutoff values — appear darkest. The points 
surrounded by white are not chaotic even for S — I. 
The cutoffs for the darkest points depend on the pa- 
rameter values: for Fig. El (a = 1 and e = 0.5), the 
same as Fig. above), we find 5'cutoff = 0.18457; Fig.lTl 
{a — 1 and e — 0.01, the same as Fig. 1341 above) has 
S-cutoff = 0.65625; and Fig. [TSl(a = 0.5 and e = 0.5, the 
same as Fig. 013) has ^cutoff ~ 0.28125. We should not 
take the digits of precision seriously, since these values 



are only accurate to within 0.05, but in all cases the spin 
cutoff values are significantly above the physically realis- 
tic range of S ^ 10~*-10~^. 

E. Retrograde orbits 

We have considered a wide variety of orbits — varying 
eccentricity and Kerr parameter for different pericenter, 
orbital inclination, and spin parameters — but so far all 
orbits have satisfied < t < 7r/2, i.e., they have all been 
prograde orbits, moving in the same direction as the cen- 
tral black hole's spin. We investigate now the case of ret- 
rograde orbits, and show that they are poor candidates 
for chaos. 

It is evident from looking at an Tp-i plot of a retro- 
grade orbit (Fig. 1161) that the pericenters are much larger 
than their prograde counterparts. For the S = 1 particle 
illustrated in Fig. 1161 the minimum empirical pericen- 
ter is larger than 6M, in contrast to prograde orbits, 
which get below 1.5 M. Furthermore, although it is clear 
from Fig. ll6f b) that the parameter space is severely dis- 
torted, the chaos is extremely weak. The largest Lya- 
punov exponent, even in this extreme S — \ case, is 
Amax = 3.5 X 10~^, two orders of magnitude smaller than 
the prograde case. Unsurprisingly, all chaos disappears 
when 5 <C 1. The smallest value of ^cutoff is 0.65265 
for the parameter values shown in Fig. we find no 
evidence of chaos below this value of S. 



F. Varying spin inclination 

So far we have always used the same values for the two 
spin components passed to the parameterization proce- 
dure (scaled by the total spin S): = = 0.2 S. We 
consider now the effect of varying these parameters, and 
also discuss the corresponding initial spin inclination an- 
gles in a fiducial rest frame. 

We begin with an initial condition that is chaotic for 
S = 1 but is otherwise arbitrary. We then vary S*'' and 
5^ and calculate the Lyapunov exponent for each con- 
figuration. The result for a — 1, ex — 0.5, rp^K — 2.3 
and LK = 20° appears in Fig. [HI When = = 0.2, 
the parameter values chosen correspond to a point from 
Fig. O There are not valid initial conditions for all 
choices of parameter; in particular, negative values of S"^ 
are often unstable or are unable to satisfy the spin con- 
straints. Nevertheless, there is a large variety of param- 
eters that do give rise to valid orbits, and many of them 
are chaotic. The strongest chaos exists for orbits that 
have small values of S^, but this appears to be mainly 
because such orbits are able to achieve small empirical 
pericenters. As Fig. ^1 clearly shows, the Lyapunov ex- 
ponent generally decreases as pericenter increases, with 
no chaotic orbits above rp^p = 2.5 M. 

While we are forced by the constraints to parameter- 
ize the equations of motion in terms of spin components. 
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(a) (b) 

FIG. 13: Spin cutoff map for o = 1 and e = 0.5. The rectangles are shaded according to the minimum value of S for which 
the corresponding initial condition is still chaotic. White rectangles indicate points that are not chaotic even when S = 1. The 
darkest points correspond to a cutoff value of S = 0.18457; below this critical value, none of the initial conditions are chaotic. 
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FIG. 14: Spin cutoff map for a = 1 and e = 0.01. The rectangles are shaded according to the minimum value of S for which 
the corresponding initial condition is stiU chaotic. White rectangles indicate points that are not chaotic even when S = 1. The 
darkest points correspond to a cutoff value of S = 0.65625; below this critical value, none of the initial conditions are chaotic. 
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FIG. 15: Spin cutoff map for a = 0.5 and e = 0.5. The rectangles are shaded according to the minimum value of S for which 
the corresponding initial condition is still chaotic. White rectangles indicate points that are not chaotic even when 5 = 1. The 
darkest points correspond to a cutoff value of S = 0.28125; below this critical value, none of the initial conditions are chaotic. 
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FIG. 16: Tp-L map of retrograde orbits (t > 90°) for a = 1, e = 0.5, 5 = 1: chaos strength as a function of pericenter and orbital 
inclination angle, (a) Requested parameters; (b) empirical parameters. The rectangles are shaded according to the Lyapunov 
exponent for the initial condition represented by each point, with darker shades of gray representing larger exponents and hence 
stronger chaos. The scaling is the same as Fig. |5] an exponent of A = 0.01 Af~^ would appear black. The chaos is weak for 
these retrograde orbits: the largest Lyapunov exponent in the plot is A = 3.5 x 10"* . 
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FIG. 17: Lyapunov exponents for varying values of initial 
spin components S"" and . The parameter values S — 1, 
Tp, K = 2.3, ek = 0.5, and ik = 20° are held fixed. The point 
with S"^ = — 0.2 appears in Fig. |5] above. The actual 
local spin inclination angles in a fiducial (ZAMO) rest frame 
appear in Table Q] Note that only one valid initial condition 
exists for negative initial S'' . 
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FIG. 19: Approximate Lyapunov exponent vs. time for a 
nonchaotic deep (rfinai = 10^ M) integration. The param- 
eter values are a = 1, 5 = 0.1, e — 0.5, Vp = 1.32 M, 
and L = 28.5°, corresponding to one of the inner orbits from 
Fig. 1281 The Lyapunov exponent appears to be zero; its time- 
evolution has the characteristic hyperbolic shape expected as 
log [re(T)]/T approaches zero for large times. A least-squares 
fit of log [T'e(T)] vs. r gives a value of A 2.8 x 10~^ M"^ . 




FIG. 18: Scatter plot of empirical pericenters r^, p vs. Lya- 
punov exponent for the spin inclinations in Fig. 1171 The Lya- 
punov exponent is primarily a function of pericenter, regard- 
less of spin inclination. 



it is more convenient to visualize the spin in a fiducial 
rest frame that is hypersurface-orthogonal to the par- 
ticle's trajectory. In Kerr spacetime, the hypersurface- 
orthogonal observers are the zero angular momentum 
observers (ZAMO), which is the same frame we use 
when calculating the Lyapunov exponents. By project- 
ing the components of the spin vector S*^ into the ZAMO 
frame in the same way as we do for the projected norm 
fSec lIVB)) . we can find the local value of the spin inclina- 
tion angle ^locai (i.e., the angle between the spin axis and 
the axis of the central black hole) . The results appear in 
Tabled It is clear that our variation of spin components 
samples a large variety of initial spin inclination angles. 



G. Deep integrations 



Since we adopted a maximum time of 10^ M for the 
Lyapunov integrations, it is reasonable to ask whether 
chaos might manifest itself on a longer timescale. This is 
certainly possible, but it appears that most initial condi- 
tions are either chaotic on a timescale of order 10^-10'' M 
or are not chaotic at all, as discussed in Sec. IIV CI An 
example appears in Fig. |21 where one initial condition is 
unambiguously chaotic, while a second located close by 
is not chaotic, even on a much longer timescale. 

To convince ourselves that slow chaos is not lurking in 
the apparently nonchaotic regions, we performed a few 
longer-time integrations. In particular, we calculated the 
Lyapunov exponents using rfinai = 10^ M for all the in- 
nermost {rp — 1.32) orbits from Fig. [2H1 (5" = 0.1) and 
Fig. I^ni {S — 10~^), which are strongly chaotic when 
S — 1 but are apparently without chaos below S = 0.1. 
The largest exponent occurs for c = 28.5°, which is 
therefore the worst-case scenario. Plots of the Lyapunov 
exponents vs. time appear in Fig. E| {S = 0.1) and 
Fig. 1201 {S ~ 10^"*); their magnitudes are on the order 
of Amax = 3 X lO^^Af^^, corresponding to an e-folding 
timescale of approximately 3.3 x 10^ M. For comparison, 
we show A vs. r for a chaotic initial condition in Fig. 1211 
it is clear that the Lyapunov exponent asymptotes to a 
nonzero value in much less than 10"" M, even for weak 
chaos. (The initial condition in Fig.[2]is the S — 0.1 or- 
bit illustrated in Figs. [TTI and [T^ whose Lyapunov expo- 
nent is actually quite small compared to analogous 5" = 1 
orbits.) The long integrations thus provide strong evi- 
dence that the disappearance of nearly all chaotic orbits 
below S* = 0.1 is a real effect. 
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TABLE 1: Local spin inclination angles Siocai in a fiducial (ZAMO) rest frame as a function of and S" . The parameter 
values S — 1, Vp^K = 2.3, Ck ~ 0.5, and lk = 20° are held fixed. An illustration of their Lyapunov exponents appears in 
Fig.CZI 
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FIG. 20: Approximate Lyapunov exponent vs. time for a 
nonchaotic deep (rflnai = lO'^ M) integration. The parameter 
values are a = 1, S = 10~^, e = 0.5, Tp = 1.32 Af, and l = 
28.5°, corresponding to one of the inner orbits from Fig. 1291 
As in Fig. 1191 the Lyapunov exponent appears to be zero. 
A least-squares fit of log [r'e(T)] vs. r gives a value of A ~ 
3.0 X 10"'^M-\ 
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FIG. 21: Approximate Lyapunov exponent vs. time for a 
chaotic deep (rflnai = lO"" M) integration. The simulation 
data is identical to that shown in Fig. 1121 the parameter values 
are a = 1, S = 0.1, ck = 0.6, r'p.x ~ 1.21 M, and lk = 
31°, corresponding to the chaotic orbit from Fig. 1321 The 
Lyapunov exponent for this chaotic orbit levels ofi' after less 
than 10^ M, in contrast to Figs. ll9l and l2UI where A continues 
to decrease even after lO'^ M. 



H. S — and chaos mimics 

Since the case of S* = corresponds exactly to 
geodesic orbits in Kerr spacetime, such systems cannot 
be chaotic — the return of the Carter constant Q and the 
loss of the spin degrees of freedom make the system fully 
integrable. Nevertheless, even some geodesic orbits can 
have a large separation of nearby initial conditions, which 
can appear to be chaotic. These chaos mimics typically 
spend many orbital periods whirling around deep in the 
strongly relativistic zone near the horizon, only occasion- 
ally zooming out to higher radii. These so-called zoom- 
whirl orbits may provide significant challenges to detec- 
tion despite their formal integrability. 

An example of how much divergence an S' = orbit 
can experience appears in Fig. |31 A picture of the cor- 



responding orbit (visualized in Boyer-Lindquist coordi- 
nates embedded in ordinary space) appears in Fig. 1221 
which makes clear the large number of low-radius 0- 
periods characteristic of zoom-whirl orbits. A second 
example of a chaos mimic appears in Figs. |^ and 1^ 
This orbit, in contrast to the previous one, does not have 
a particularly small pericenter, but its high inclination 
angle and zoom-whirliness allow it to mimic chaotic or- 
bits. 

The chaos mimics can exhibit large growth of the ini- 
tial deviation vector, approximately a factor of 10^-10*, 
after 10^ M. The principal means for detecting them 
is by requiring several high-separation points in a row 
(on a time-T « 100 M basis) as mentioned in Sec. IIV CI 
the mimics have oscillations with high amplitudes due 
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(a) (b) (c) 

FIG. 22: The orbit of a chaos mimic: a non-spinning {S = 0) test particle in maximal (a = 1) Kerr spacetime, plotted in 
Boyer-Lindquist coordinates. A Lyapunov plot of these initial conditions appears in Fig. |21 from Sec. liV CI (a) The orbit 
embedded in three-dimensional space, treating the Boyer-Lindquist coordinates as ordinary spherical polar coordinates; (b) 
y = r sin 9 sin vs. a; = r sin 6 cos (f>; (c) z vs. p = ^ ^'ip-. The inclination angle is t = 31° , while the pericenter is Tp = 1.1 M 
(just 0.1 AI above the horizon at r_r/ = 1 M). 




FIG. 23: Another chaos mimic: a non-spinning {S = 0) test particle in maximal (a = 1) Kerr spacetime, plotted in Boyer- 
Lindquist coordinates. A Lyapunov plot of these initial conditions appears in Fig. 1241 (a) The orbit embedded in three- 
dimensional space, treating the Boyer-Lindquist coordinates as ordinary spherical polar coordinates; (b) y = r sin 6 sin (j) vs. 
X = r sin ^ cos 0; (c) z ys. p — y'x^ + y^. The inclination angle is t = 88.5°, while the pericenter is rp = 4.4 Af. 



to their zoom-whirliness, but they do not represent true 
saturations of the separation vector. 

VI. CONCLUSIONS 

The Papapetrou equations, which model a spinning 
test particle, exhibit chaotic solutions in Kerr spacetime 
for a wide range of parameters. In terms of the mass M of 
the central (Kerr) black hole, the largest Lyapunov expo- 
nents are of order A„iax = 0.01 M~^, which represents an 
exponential divergence of trajectories on a timescale of 
T\ = 1/X = 100 M. Furthermore, there are many chaotic 



orbits with exponents in the range 10~'^-10~^ M~^. De- 
spite the large number of chaotic orbits, we find that 
values of A corresponding to unambiguous chaos occur 
exclusively when the spin parameter S is not small com- 
pared to unity. In particular, we find virtually no chaos 
for spin values below S = 0.1, and no evidence of any 
chaos for spins below S = 10^'*. 

The strongest determinant of chaotic behavior, apart 
from the spin parameter, is the pericenter of the orbit in 
question. The most highly chaotic orbits are those that 
reach pericenters near the horizon of the black hole. This 
is due to the high spacetime curvature in these regions, 
which maximizes the size of the coupling of the spin to 
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FIG. 24: A chaos mimic: the natural logarithm of the prin- 
cipal ellipsoid axis (log [7'e(T)]) vs. r for an S — orbit. The 
size of the initial deviation vector is eo = 10~*. The value of 
log [re(T)] briefly rises up to the saturation level at log(0.9/eo) 
[so that e(T) = 0.9], but the orbit is not chaotic since its spin 
satisfies S = 0. A plot of the corresponding orbit appears in 
Fig. 113 



APPENDIX: SPIN VECTOR FORMULATION 

We summarize here the formulation of the Papapetrou 
equations in terms of the spin 1-form S*^ (often referred 
to loosely as the "spin vector"), as mentioned in Sec. Ill Al 
In this paper we consider a spinning particle of rest 
mass fj, orbiting a central Kerr black hole of mass M, 
and it is convenient to measure all times and lengths in 
terms of M and all momenta in terms of /i. In these nor- 
malized units, the equations of motion in terms of the 
spin 1-form are 

= v'^ 

dr 

ViTP,, = -Rl^^v^p^Sp (A.l) 

where 



the Riemann curvature tensor [Eq. (|2.H . When the Kerr 
parameter a is small, so that the Kerr metric differs only 
slightly from the Schwarzschild metric, the pericenters 
are much higher than in the extreme Kerr (a = 1) case. 
Chaos in the Papapetrou system is therefore weak when 
a is small. The prevalence of chaos is also dependent 
on orbital eccentricity. Near-circular (e = 0.01) orbits 
have many fewer regions of chaotic orbits than those with 
higher eccentricities (e = 0.5 or e = 0.6). This seems 
due primarily to the lower pericenters accessible to high- 
eccentricity orbits. 

The dependence of the Lyapunov exponents on S is our 
most important result: in all cases considered, physically 
realistic values of S (satisfying S 1) are not chaotic. 
We have shown conclusively that the Papapetrou equa- 
tions admit many solutions that arc formally chaotic, but 
without exception such chaotic solutions occur only for 
relatively large values of S. Below the upper limit for 
physically realistic spins {S ^ 10~^), we find no evidence 
of chaotic solutions. As a practical matter, this means 
that chaos will not manifest itself in the gravitational 
radiation from extreme mass-ratio binary inspirals. 

Acknowledgments 

I would like to thank Scott Hughes and Teviet 
Creighton for providing notes on parameterizing Kerr 
geodesies in terms of orbital parameters. I also thank 
Sterl Phinney for his careful reading of the manuscript 
and perceptive comments. This work was supported in 
part by NASA grant NAG5-10707. 



= iRV"^"- (A-2) 

The supplementary condition Eq. H2.5|l allows for an 
explicit solution for the 4- velocity in terms of p^: 

vf" ^ Nip" + W''), (A.3) 

where 

^ -*R*^^'^P■ys^Pf}S^ (A.4) 

and 

The normalization constant N is fixed by the constraint 
v^vi" = -1. 

The spin 1-form satisfies two orthogonality constraints: 
P^Sf, = (A.6) 

and 

v^'S^ = 0. (A.7) 

These two constraints are equivalent as long as is given 
by Eq. (ESJ: = yf^Sf^pf'S^ + w'^Sf, =p''S^, since by 
definition of w^^ [Eq. (|A.4|l ] the second term involves the 
contraction of a symmetric tensor with an antisymmetric 
tensor and therefore vanishes. We enforce Eq. (jA.6p in 
our parameterization scheme, and we use Eq. lfO)l in the 
equations of motion, so Eq. IIA.7|) is then automatically 
satisfied. 
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FIG. 25: Vp-L map for S = 0.9, a = 1, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. 1^1 Chaotic orbits are widespread. 
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FIG. 26: Tp-t map for S = 0.5, a = 1, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. El Because of the extremely low pericenters accessible at this value of 5* 
(which are excluded when S = 1), the chaos for S — 0.5 is the strongest we find. The largest Lyapunov exponent is just over 
A = 0.01 A/~\ 
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FIG. 27: rp-L map for S = 0.2, a — 1, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. ^ Only a few orbits are chaotic. 
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FIG. 28: Vp-L map for S — 0.1, a — 1, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is 
scaled to the same maximum Lyapunov exponent as Fig. |3| All chaos is gone, although the parameter space is still somewhat 
compressed. 
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FIG. 29: 



-L map for 5* = 10 ^, a = 1, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is 



scaled to the same maximum Lyapunov exponent as Fig. ^ There are no chaotic orbits. The empirical parameter values in (b) 
are indistinguishable from the requested values except for initial conditions that specify values of t corresponding to unstable 
orbits (as discussed in Sec. 1111 B 3ll . 
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FIG. 30: rp-L map for S = 1, a = 1, and e = 0.6. (a) Requested parameters; (b) empirical parameters. The shading is 
scaled to the same maximum Lyapunov exponent as Fig. El Chaotic orbits are widespread. The largest Lyapunov exponent is 
A = 5.5 X 10"^ M-^ 
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FIG. 31: rp-L map for S = 0.5, a = 1, and e = 0.6. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. |5| As in the case of Fig. 1261 the low accessible pericenters give rise to strong 
chaos. The largest Lyapunov exponent is A = 9.2 x 10~"^M~^. 
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FIG. 32: rp-L map for S — 0.1, a — 1, and e = 0.6. (a) Requested parameters; (b) empirical parameters. The shading is 
scaled to the same maximum Lyapunov exponent as Fig. ^ There is only one chaotic initial condition (which is in fact the 
only 5 = 0.1 chaos we find), but the chaos is real, as discussed in Sec. IV Bl and illustrated in Fig. 1121 The Lyapunov exponent 
is A = 1.0 X 10"^ M-^ 
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FIG. 33: Tp-t map for 5 = 10~*, a = 1, and e — 0.6. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig.|3] All chaos has disappeared. As in Fig. 1291 the empirical parameter values 
in (b) are indistinguishable from the requested values except for initial conditions that specify unstable orbits (Sec. 1111 B 3ll . 
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FIG. 34: 



map for S = 1, a — 1, and e = 0.01. (a) Requested parameters; (b) empirical parameters. The shading is scaled 



to the same maximum Lyapunov exponent as Fig. |5] There is some relatively weak chaos. The largest Lyapunov exponent is 
A = 2.7 X 10"^ M" 
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FIG. 35: Vp-L map for S — 0.5, a = 1, and e = 0.01. (a) Requested parameters; (b) empirical parameters. The shading is 
scaled to the same maximum Lyapunov exponent as Fig.|5| There are a few regions of weak chaos. 
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FIG. 36: Vp-L map for S = 0.1, a = 1, and e = 0.01. (a) Requested parameters; (b) empirical parameters. The shading is 
scaled to the same maximum Lyapunov exponent as Fig.|5| There are apparently no chaotic initial conditions. 
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FIG. 37: 



-L map for S — 10 



1, and e — 0.01. (a) Requested parameters; (b) empirical parameters. The shading is 



scaled to the same maximum Lyapunov exponent as Fig. |5| The chaos had disappeared. 
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FIG. 38: Vp-L map for S = 1, a — 0.9, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. 1^1 Chaotic orbits are widespread. 
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FIG. 39: Tp-i map for S = 1, a = 0.7, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The sliading is scaled 
to the same maximum Lyapunov exponent as Fig. 1^1 There is still substantial chaos. 
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FIG. 40: Tp-i map for S = 1, a = 0.5, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. El The chaos is largely confined to low pericenter orbits. 
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FIG. 41: Tp-L map for S = 1, a = 0.4, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled 
to the same maximum Lyapunov exponent as Fig. |2| Only a handful of initial conditions are chaotic. 
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FIG. 42: rp-L map for S = 1, a — 0.2, and e = 0.5. (a) Requested parameters; (b) empirical parameters. The sliading is scaled 
to tlie same maximum Lyapunov exponent as Fig. 1^1 No initial conditions are chaotic. 
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FIG. 43; rp-i map for S = 1, a = 0, and e — 0.5. (a) Requested parameters; (b) empirical parameters. The shading is scaled to 
the same maximum Lyapunov exponent as Fig. El No initial conditions are chaotic. Note that every column of (a) is identical. 
This is a result of the spherical symmetry of the a = (Schwarzschild) metric: all inclination angles are equivalent. As seen 
in (b), this symmetry is broken by the spin of the test particle. 



